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ABSTRACT 

Context. Thermal emission by stochastically heated dust grains (SHGs) plays an important role in the radiative transfer (RT) problem 
for a dusty medium. It is therefore essential to verify that RT codes properly calculate the dust emission before studying the effects of 
spatial distribution and other model parameters on the simulated observables. 

Aims. We define an appropriate problem for benchmarking dust emissivity calculations in the context of RT simulations, specifically 
including the emission from SHGs. Our aim is to provide a self-contained guide for implementors of such functionality, and to 
offer insights in the effects of the various approximations and heuristics implemented by the participating codes to accelerate the 
calculations. 

Methods. The benchmark problem definition includes the optical and calorimetric material properties, and the grain size distributions, 
for a typical astronomical dust mixture with silicate, graphite and PAH components; a series of analytically defined radiation fields 
to which the dust population is to be exposed; and instructions for the desired output. We process this problem using six RT codes 
participating in this benchmark effort, and compare the results to a reference solution computed with the publicly available dust 
emission code DustEM. 

Results. The participating codes implement different heuristics to keep the calculation time at an acceptable level. We study the effects 
of these mechanisms on the calculated solutions, and report on the level of (dis)agreement between the participating codes. For all but 
the most extreme input fields, we find agreement within 10% across the important wavelength range 3 gm < A < 1000/mi. 
Conclusions. We conclude that the relevant modules in RT codes can and do produce fairly consistent results for the emissivity 
spectra of SHGs. This work can serve as a reference for implementors of dust RT codes, and it will pave the way for a more extensive 
benchmark effort focusing on the RT aspects of the various codes. 

Key words. Radiation mechanisms: thermal - dust, extinction - Infrared: ISM - Radiative transfer - Methods: numerical 


1. Introduction 

Dust substantially affects the radiation emerging from many as- 
trophysical systems. To study the three-dimensional structure of 
these systems, it is often useful to numerically simulate the trans¬ 
port of radiation through a model that includes a dusty medium 
with appropriate characteristics. Many authors have described 
radiative transfer (RT) codes designed to tackle this problem; 
for an overview see for example Whitney (2011) and Steinacker 
et al. (2013). In multi-wavelength studies, thermal emission by 
the dust plays an important role. While larger dust grains can of¬ 
ten be assumed to be in local thermal equilibrium (LTE) with the 


surrounding radiation field, this assumption does not typically 
holds for very small grains (VSGs) or for polycyclic aromatic 
hydrocarbon molecules (PAHs). The absorption of a single op¬ 
tical photon substantially boosts the internal energy of such a 
small collection of atoms, causing its emission spectrum to vary 
over time (see, e.g., Sellgren 1984; Dwek 1986; Boulanger & 
Perault 1988; Helou et al. 2000; Draine et al. 2007; Smith et al. 
2007). We use the term stochastically heated grains (SHGs) to 
collectively indicate dust grains and PAH molecules that can not 
be assumed to be in LTE with the radiation field. 
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Since SHGs spend a significant amount of time at much 
higher energy levels than if they would be in LTE, they emit 
at shorter wavelengths, which can have an important effect on 
the observed spectmm. It is therefore essential to verify that RT 
codes properly calculate this emission before studying the effects 
of spatial distribution and other model parameters on the simu¬ 
lated observables. In this paper we present a benchmark prob¬ 
lem for this purpose, including a dust model, a series of input 
radiation fields, and a reference dust emission spectrum for each 
input field. The dust model described here has been designed for 
use in the TRUST benchmarks 1 (Transport of Radiation through 
a dUSTy medium), which test the actual RT aspect of various 
codes. 

A typical 3D RT simulation calculates the dust emission 
spectra for millions of dust cells (or at least for many thousands 
of library items that are representative of the cells). When cal¬ 
culating the emission spectrum for the dust grain population in a 
particular cell, the first task is to determine the temperature prob¬ 
ability distribution of the grains, given the grain sizes and chemi¬ 
cal compositions, and given the radiation field in the cell. This is 
also the computationally most demanding part of the calculation. 
In the current epoch, performing a full treatment of vibrational 
quantum modes for the dust grains in each cell is computation¬ 
ally prohibitive. In practice, RT codes use the continuous cooling 
assumption, which was shown to provide a good approximation 
in areas relevant for RT through the interstellar medium (Draine 
& Li 2001). 

The reference solutions in this paper were produced with the 
public version of the DustEM code (Compiegne et al. 2011). 
DustEM determines the grain temperature distribution by iter¬ 
atively solving an integral equation (Desert et al. 1986). We 
compare the reference solutions with the dust emission spectra 
calculated by six distinct RT codes. These codes determine the 
grain temperature distribution by solving a set of linear equations 
(Guhathakurta & Draine 1989). While this method is inherently 
faster (Guhathakurta & Draine 1989), it still becomes very ex¬ 
pensive when the grains are in LTE with the radiation field. To 
avoid this problem, care must be taken to properly transition the 
calculation from the stochastic to the equilibrium regime. 

More generally, the need for fast dust emission calculations 
has prompted the authors of RT simulation codes to implement 
various acceleration techniques, discretization choices, approx¬ 
imations, and heuristics. We study the effects of these mech¬ 
anisms on the calculated solutions, and quantify the level of 
(dis)agreement between the participating codes, with the objec¬ 
tive to help inform the interpretation of RT simulation results that 
include SHG dust emission calculations of the type presented 
here. 

The information provided in this paper and in the accompa¬ 
nying data files 2 is self-contained. Readers can implement the 
code to calculate the SHG emission, set up the benchmark tests, 
and verify the results, solely based on the information provided 
here, without referring to other sources. 

In Sect. 2 and Sect. 3 we define the dust model and the in¬ 
put radiation fields used in the benchmark problem. Section 4 
presents the linear equation method used by the RT codes repre¬ 
sented in this paper to calculate the SHG emission spectra, and 
Sect. 5 briefly introduces each of these codes. In Sect. 6 we com¬ 
pare the results produced by the RT codes with the reference so¬ 
lution, and we discuss the differences between the methods used 


1 TRUST benchmarks: http://ipag.osug.fr/RT13/RTTRUST/ 

2 SHG benchmark data: http: //www. shg. ugent. be 


by the various codes and how they influence the results. Finally, 
in Sect. 7 we summarize and conclude. 

For ease of reference. Table 1 lists the symbols used in this 
paper for various physical quantities, with the corresponding SI 
units. 


2. Dust model 

The exact choice of dust grain model (e.g., Weingartner & 
Draine 2001a; Zubko et al. 2004; Jones et al. 2013) is not crit¬ 
ical for benchmark purposes, as long as all codes employ the 
same model. Specifically, our choices do not imply a preference 
for or an endorsement of a particular model. Still, to properly 
evaluate the results of our calculations and the effects of ap¬ 
proximations in the context of real-world RT simulations of as- 
trophysical systems, it is important to use grain compositions 
and size distributions that are in agreement with observational 
constraints, rather than defining arbitrary synthetic grain prop¬ 
erties. With this in mind, we have elected to utilize the simple 
BARE-GR-S model of Zubko et al. (2004). This model uses a 
mixture of spherical, single composition (BARE) graphitic (GR) 
and silicate (S) dust grains. The graphitic grains include both 
graphite grains as well as astronomical PAH molecules. The rel¬ 
ative populations of the three components and their size distribu¬ 
tion has been optimized to reproduce observations characteristic 
of the diffuse Milky Way interstellar medium with Ry = 3.1. The 
observational constraints on the model include the shape of the 
wavelength dependent extinction, the infrared emission, and the 
abundance of elements locked up in the solid phase. 

Assuming that the scattering process is elastic, so that the 
internal energy of the interacting dust grain remains unaffected, 
scattering is irrelevant for the benchmark problem presented in 
this paper. Thus the subsequent discussion does not focus on the 
scattering properties of the dust model. 

2.1. Optical grain properties 

For both the graphite and silicate components, optical proper¬ 
ties were computed directly from the dielectric functions (re¬ 
fractive indices) of the bulk material using a Mie code originally 
provided by Viktor Zubko with some small modifications. The 
optical properties include absorption efficiencies Qf ls (a, A) as a 
function of composition k, grain size a, and wavelength A. 

While the dielectric functions are in general functions of both 
temperature and size (Draine & Lee 1984), to minimize free pa¬ 
rameters and in keeping with essentially all astronomical appli¬ 
cations of grain properties, we adopt a single set of dielectric 
functions at a specific temperature and size for each component. 
The graphitic indices were taken from Draine (2003) for 0.1 pm 
grains at 20 K. Calculations were carried out for both parallel 
and perpendicular orientations relative to the basal plane and 
combined using the 1/3-2/3 approximation (Draine 1988). For 
the silicate grains, the dielectric functions for smoothed astro¬ 
nomical silicates (0.1 pm) (Laor & Draine 1993; Weingartner 
& Draine 2001a) were used in the Mie calculations. With these 
sets of dielectric functions, assuming sperical grains, the effi¬ 
ciencies can be calculated for each component as a function of 
size and wavelength. For each component, the Mie calculation 
was carried out for 121 logarithmically spaced sizes between 
0.00035 pm and 100 pm. For each size, optical properties were 
computed at 1201 wavelength points logarithmically spaced be¬ 
tween 0.001 pm and 10000 pm. The efficiencies thus derived can 
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Table 1. Symbols used in this paper for various physical quantities, with the corresponding SI units. 


Symbol 

Units 

Description 

A 

m 

Wavelength 

s 

m 

Distance along a path 

V 

m 3 

Volume 

M 

kg 

Dust mass 

T 

K 

Temperature 

T 

s 

Interaction timescale 

p 

kgm~ 3 

Dust mass density 

£j_abs,sca,ext 

m 2 

Cross section (absorption, scattering or extinction) 

A4t 

H 

Number of hydrogen atoms 

«H = Wh/V 

Hnr 3 

Hydrogen atom number density 

P = M/N h 

kgH 1 

Dust mass per hydrogen atom 

g = ct/Nh 

m 2 H 1 

Cross section per hydrogen atom (absorption, scattering or extinction) 

k = cr/M = g/ t u 

nrkg- 1 

Mass coefficient (absorption, scattering or extinction) 

a 

m 

Dust grain size 

n d 

i 

Number of dust grains 

«d = Nd/V 

m- 3 

Dust grain number density 

«(a) = (i?)/*h 

m" 1 H- 1 

Dust grain size distribution per hydrogen atom 

^abs,sca,ext 

1 

Efficiency (absorption, scattering or extinction) 

pbulk 

kgm~ 3 

Bulk mass density of grain material 

C 

JK-'kg- 1 

Specific heat capacity of grain material 

h 

Jkg- 1 

Specific enthalpy (internal energy per unit mass) 

H 

J 

Enthalpy (internal energy) of a dust grain 

J 

W m -3 sr 1 

Mean spectral radiance (intensity of radiation field) 

B 

W m~ 3 sr _1 

Black-body spectral radiance (Planck’s law) 

U 

1 

Radiation field strength relative to solar neighborhood 

j 

Wnr 1 sr 1 

Spectral dust emission 

s = j/Mt 

Wm-'sr'H - 1 

Spectral dust emission per hydrogen atom 


be cast as either cross sections or mass coefficients: 


2.2. Grain size distributions 


abs.sca/ i\ 2/^abs,sca 

cr ( a , A) = na Q 


abs.sca/ v\ 

k, ( a , A) = 


32 ; 


k 

.abs.sca 


(a, /I) 
{a, /l) 


4 ap k 


(1) 

( 2 ) 


As there are no refractive indices for the class of materials 
commonly referred to as astronomical PAH, the optical proper¬ 
ties for these materials are simply constructed so as to reproduce 
the available observations. Therefore, the PAH cross sections 
were computed following Draine & Li (2007). The utilization 
of the Draine & Li (2007) formulation for the cross sections rep¬ 
resents a deviation from a pure Zubko et al. (2004) BARE-GR-S 
model as Zubko et al. (2004) utilized the Li & Draine (2001) 
PAH cross sections. 

In Draine & Li (2007), the PAH cross sections were up¬ 
dated to reflect knowledge gained regarding the shape of the 
PAH emission spectrum from the wealth of Spitzer observations 
available. Also of note, the PAH component in the model defined 
here consist of pure neutral PAHs; no attempt was made to ac¬ 
count for the varying ionization fraction of a PAH molecule as a 
function of effective PAH size and radiation field. The PAH op¬ 
tical properties were computed on the same wavelength grid as 
the graphite and silicate components of the model. However, the 
size grid was of course altered; properties were computed for 28 
logarithmically spaced sizes between 0.00035 p m and 0.006 pm. 


The relative contribution of each material in a dust model is 
set by the grain size distribution for that particular material. 
The shape of the size distributions in the BARE-GR-S model 
of Zubko et al. (2004) matches observations of the interstel¬ 
lar medium in the solar neighborhood, including the amount of 
refractory material available, as well as the wavelength depen¬ 
dence of the Milky Way diffuse (Ry = 3.1) extinction and emis¬ 
sion spectrum. Zubko et al. (2004) provide convenient analytic 
functional approximations to the size distributions for each com¬ 
ponent of the model, which take the form 


log O(a) = c 0 + b 0 log 



■b i 


log 

f ? 3 


/ a 
\0i 
Cl — Cl 3 


bi 


log - 
\a 2 


/um 


■ b 4 


a - a 4 


yum 


(3) 


where a is the grain size and a,, /?,, c,, m l represent constant pa¬ 
rameters. In the interest of clarity, the notation in Eq. 3 omits a 
subscript k indicating the type of material being referenced. For 
the BARE-GR-S model, we require three sets of parameters, one 
for each of the materials comprising the dust model. The appro¬ 
priate values can be found in Table 7 of Zubko et al. (2004), and 
are reproduced here in Table 2. 
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Table 2. Parameters for the analytical approximation to the BARE-GR-S size distribution defined in Eqs. 3 through 5 and the bulk densities for the 
three components in our dust model. With these parameter values, the grain size a substituted in the equations must be expressed in pm and the 
resulting size distribution £l k (a) is expressed in number of dust grains per hydrogen atom and per pm. 


Parameter 

Units 

PAH 

Graphite 

Silicate 

a min 

fim 

0.00035 

0.00035 

0.00035 

a max 

li m 

0.005 

0.33 

0.37 

A 

/tm- 1 H- 1 

2.227433X10- 7 

1.905816xl0- 7 

1.471288X10- 7 

co 

1 

-8.02895 

-9.86000 

-8.47091 

bo 

1 

-3.45764 

-5.02082 

-3.68708 

b\ 

1 

1.18396xl0 3 

5.81215xl0- 3 

2.37316xl0- 5 

a\ 

jjm 

1 

0.415861 

7.64943xl0- 3 

m\ 

1 

-8.20551 

4.63229 

22.5489 

bi 

1 

0 

0 

0 

as 

fjm 

1 

1 

1 

m 2 

1 

1 

1 

1 

b 3 

1 

l.OxlO 24 

1.12502xl0 3 

2.96128xl0 3 

03 

jjm 

-5.29496xl0- 3 

0.160344 

0.480229 

m3 

1 

12.0146 

3.69897 

12.1717 

b A 

1 

0 

1.12602xl0 3 

0 

04 

/im 

0 

0.160501 

0 

m 4 

1 

1 

3.69967 

1 

pbulk 

kg nG 3 

2240 

2240 

3500 


By construction, Eq. 3 has been normalized such that 

r r 

O k (a)dci = 1 (4) 

Jaf n 

where a™ n and a™ ax are the minimum and maximum sizes over 
which the size distribution for component k is defined. The com¬ 
plete size distribution function is then given by 

D. k (a) = A k O k (a ) (5) 

where A k is the overall normalization of component k, as listed 
in Table 2. 

With this definition of the grain size distribution, the total 
number of dust grains per hydrogen atom and the total dust mass 


per hydrogen atom are given by: 


Afo/M i = Yj f ^k(a)da 

k Ja T 

(6) 

p = M/N h = ^ f ^7ra 3 p bulk Q4a) da. 

k Ja T° 3 

(7) 


We also provide the grain size distributions defined by Eqs. 3 
and 5 and Table 2 in tabulated form (see Sect. 2.4), on a size grid 
with 24, 62, and 63 samples for the PAH, graphite, and silicate 
component, respectively. 

2.3. Calorimetric grain properties 


The internal energy of a dust grain of composition k and its tem¬ 
perature are related via the heat capacity of the grain material: 


H k (a,T) = 

^7ra 3 p bu %.(7’) 

r T 

(8) 

h k (T ) = 

c k (T')dT' 

Jo 

(9) 


where H(a, T ) is the internal energy (enthalpy) of a grain with 
size a at temperature T; h(T ) is the specific enthalpy (per unit 
mass); c{T) is the specific heat capacity; and p bulk is the bulk 
density of the grain material. 

We elected to use the heat capacity functions proposed in 
Draine & Li (2001) and Li & Draine (2001). To avoid subtle dif¬ 
ferences between implementations, we provide this information 
in tabulated form (see Sect. 2.4) rather than expecting each code 
to implement the equations. One table describes the graphitic 
components (PAH molecules and graphite grains) and another 
table describes the silicate component. Each table lists the spe¬ 
cific enthalpy and the specific heat capacity at 1000 logarithmi¬ 
cally spaced temperature points ranging from 1 K to 2500 K. 

Finally, the bulk densities for the three components in our 
dust model are listed in Table 2. 

2.4. Data files 

The data files defining the dust properties described above can 
be downloaded from the web site indicated in footnote 2. They 
are contained in the DustModel directory, which is organized in 
the following subdirectories: 

- Grainlnputs: optical and calorimetric properties for each 
of the three grain types k. Optical properties include the ab¬ 
sorption efficiencies <2 abs (cz,/l) tabulated on a grid of 1201 
wavelengths A and 28 (for PAH) or 121 (for graphite and sil¬ 
icate) grain sizes a. Calorimetric properties include the spe¬ 
cific enthalpy h(T) and the specific heat capacity c(T) of the 
grain material, tabulated on a grid of 1000 temperatures T. 
The calorimetry data for graphitic grains should be used for 
PAH molecules as well. 

- Sizelnputs: tabulated grain size distributions fl k (a) for 
each grain type k, on a size grid with 24, 62, and 63 samples 
for the PAH, graphite, and silicate component, respectively. 
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Implementations may choose to compute the size distribu¬ 
tion from the functional form defined in Sect. 2.2, or to load 
the tabulated data. 

- EffectiveGrain: size-integrated values of the optical 
properties. This information is not needed for calculating 
dust emission. 

- ScatMatrix: scattering matrix elements. This information 
is not needed for calculating dust emission. 


3. Radiation fields 

3. 1. Basic definitions 

The spectral radiation for a black body in thermal equilibrium at 
temperature T in function of wavelength A is given by the Planck 
function 


B(A, T) 


2hc 2 1 


A 5 


exp( 


he 

AkT 


)-i 


( 10 ) 


where li denotes the Planck constant, c the speed of light in vac¬ 
uum, and k the Boltzmann constant. 

We define the solar neighborhood interstellar radiation field 
(ISRF) given in table A3 of Mathis et al. (1983) through the 
following functional form inspired by (but not identical to 3 ) Eq. 
31 in Weingartner & Draine (2001b): 


7 Mat (d) = 


0 

A < 0.0912 pm 

3069 W/m 3 /sr x (A/p m) 3 4172 

0.0912/un < A < 0.110/mi 

1.627 W/m 3 /sr 

0.110/un < A < 0.134/mi 
0.0566 W/m 3 /sr x (A/p m)" 1 - 6678 

0.134/un < A < 0.250 /on 
10~ 14 B(/1,7500K) + 1CT 13 B(/l, 4000K) 

+ 4 x 1CT 13 B(A, 3000 K) 

0.250/rm < A 


( 11 ) 


Note that the recipes in the other papers prescribe the total radi¬ 
ation field 4^7 Mat , whereas we prescribe the mean radiation field 
7 Mat . Based on this reference field, we define the strength U of 
an arbitrary radiation field J(A) as 


u= f J(A)dA 1 

/-»oo 

y Mat (d)dd 

(12) 

Jo 1 

Jo 



3.2. Benchmark input fields 

In the benchmark described in this paper, the dust grains are ex¬ 
posed to two sets of distinct radiation fields. The first set consists 
of eleven scaled versions of the Mathis ISRF, ranging from weak 
to strong , defined as 

i SHG ’'(T) = Ui x 7 Mat (/l) with U, = 1(T 4 ,1(T 3 ,..., 10 5 ,10 6 (13) 

3 The Weingartner & Draine (2001b) equation is formulated in func¬ 
tion of frequency rather than wavelength, and the dilution factor for the 
4000 K black body listed in their Table 1 is not adjusted to the value 
specified in Sect. 2.1 of Mathis et al. (1983). 


The second set consists of the following six diluted black body 
fields with varying temperatures, ranging from soft to hard: 

8.28 x 1(T 12 B(A, 3000 K) 

2.23 x 1CD 13 B(/l, 6000K) 

2.99 x \Qr u B(A, 9000 K) 

7.23 x 1(T 15 B(A, 12000 K) 

2.36 x 1(T 15 B(A, 15000 K) 

,9.42 x 1(T 16 £(4,18000 K) 

The dilution factors were chosen so that the far-infrared peak of 
the dust emissivity is at the same level for all fields in this set 
(for ease of visualization), and so that all fields in the set have a 
strength of 1 ^ U < 10. 

3.3. Calculation and wavelength grid 

Using the dust model described in section 2, the codes partici¬ 
pating in this benchmark calculate the spectral dust emissivity 
s(A) for each of the input radiation fields specified in section 3, 
taking into account stochastic heating of small grains. The radia¬ 
tion emitted by the dust itself is ignored with respect to the input 
field, i.e. it is not the intention to calculate a self-consistent radia¬ 
tion field. The calculations are performed, and the results written 
down, using the wavelength grid on which the optical properties 
have been tabulated. This is a logarithmic grid with 1201 points 
in the range 0.001 pm < A < 10000 pm. 


J smj (A ) = < 


4. Dust emission 

4. 1. Emission from a dust mixture 


The thermal emission of a dust grain depends nonlinearly on the 
grain size a. even in LTE conditions. It is therefore impossible 
to calculate the emission for a dust mixture with varying grain 
sizes from effective grain properties that would somehow rep¬ 
resent the whole mixture (Steinacker et al. 2013; Wolf 2003). 
Instead, we define a grid of grain size bins b for each dust model 
component k, and we choose an average, representative grain for 
each bin. We then proceed to calculate the emission as if each bin 
would contain only representative grains. For a sufficiently large 
number of bins, this procedure converges to the proper result. 

A simple approach is to represent each bin by a grain size at 
the arithmetic or geometric center of the bin. In a somewhat more 
sophisticated approach, the absorption cross section per hydro¬ 
gen atom representative for a particular bin can be calculated by 
an integration over the size distribution: 




ttec Qf’ s (a, d) Clk(a) do 


(15) 


where [a™f, a™ x ] specifies the size range of bin b for dust model 
component k. 

The representative mass of a dust grain in a particular bin can 
similarly be obtained from 

r a kf 

I ' Sl k (a)da (16) 

so that the enthalpy of a representative dust grain at temperature 
T is given by 


Mk,b 


f^k.b 

aft” 


A jr 

Pit ulk Qk(a) do 


H k , h (T) = M k b h k (T) 


(17) 
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where h^T) is the specific enthalpy of the grain material of dust 
component k at temperature T. 

The emissivity per hydrogen atom from a dust mixture with 
grain type components k and grain size bins b exposed to a radia¬ 
tion field J(A), called the input field, can be expressed in function 
of the representative grain properties as 

r*co 

eU) = V P k , hJ (T) B(A, T) d T (18) 

k* ' Jo 


where B(A, T) is the Planck function defined in Eq. 10, and 
Pk,bj(T) is the probability of finding the representative grain of 
bin k, b at temperature T. 

The emission originating from a dust mixture with specified 
total mass M can then be written as 

j(A) - — s(A) (19) 

B 

with p given by Eq. 7. When combining the emission from var¬ 
ious dust mixes, it is useful to recall that it is physically mean¬ 
ingful to add cross sections and masses, while mass coefficients, 
in general, cannot be added meaningfully: 


^ ft ft 

K\ + K2 - - + - 

Pi P2 


* 


ft + ft 
Pi +P2’ 


( 20 ) 


The challenge is thus to compute the probability distribution 
of grain temperatures, Pk.bj(T), which depends on the input ra¬ 
diation field in addition to the grain properties. See, for example. 
Fig. 4 of Draine & Li (2007) for an illustration of various tem¬ 
perature distribution curves. 

In this discussion, we characterize P as a function of grain 
temperature. The temperature of a grain and its internal energy 
are related through Eq. 8, so we could equivalently characterize 
P as a function of internal grain energy. 


4.2. Equilibrium heating dust emission 

When the representative grain in bin k, b is in LTE with the sur¬ 
rounding radiation field J(A), the temperature probability distri¬ 
bution Pk.b.j(T) becomes a delta function at the grain equilibrium 
temperature 

Pk.b.j(T) = SOT ~ THj) (21) 

and the equilibrium temperature can be determined via the en¬ 
ergy balance equation 


a specific radiation field, dropping the indices k, b and J from the 
notation. 

We select an appropriate temperature grid with N bins 7), 
i = 0,... ,N - 1 (see Sect. 4.1). Our goal is to determine the 
probabilities P, = P(P,) for a grain to reside in temperature bin 
i. 

We define a transition matrix A fj describing the probability 
per unit time for a grain to transfer from initial temperature bin i 
to final temperature bin /. The transition matrix elements in the 
case of heating (/ > i) are given by 

, he AH r 

A fj - 4 7T <? bs (A fi ) J(A fi ) --^—3 (23) 

[H(T f ) - H(Ti)\ 

where H(Tf) and H(Tj) are the enthalpies of a dust grain in the 
final and initial temperature bins, A Hf = //(T alax ) - //(T" lm ) is 
the enthalpy width of the final temperature bin, and Aft is the 
transition wavelength which can be obtained from 

he 

A fi = -. (24) 

; H(T f ) - H(Tj) 

We assume that a dust grain cools by radiating photons with an 
energy that is very small compared to the internal energy of the 
grain. With this continuous cooling approximation, cooling tran¬ 
sitions occur only to the next lower level, so that A fj - 0 for 
/ < i — 1 and 


A 


/-i,i 


47T 


H(Ti) - urn ,) 



<T abs (/l) B(A, Tj) d A. 


(25) 


The diagonal matrix elements are defined as A ti = - 2 /#i A fj, 
however there is no need to explicitly calculate these values as 
they are not used in the final procedure. 

Assuming a steady state situation, the probabilities P, can be 
obtained from the transition matrix by solving the set of N linear 
equations 


N -1 

Yj A fjPi = 0 / = 0, ...,N - 1 (26) 

;=o 

along with the normalization condition 


v-1 



(27) 


/"»CO /-»00 

gl b f(A)J(A)dA= ^(A)B(A,T eq bJ )dA. (22) 

Jo Jo 

4.3. Stochastic heating dust emission 

When a single photon absorption may significantly change the 
internal energy of a representative grain, the grain is not in LTE 
with the surrounding radiation field J(A). The grain is stochas¬ 
tically heated, and its state can no longer be characterized by a 
single temperature. In that case, we need to solve for the tem¬ 
perature probability distribution Pk,b,j(T ) to calculate the grain 
emission. The six RT codes benchmarked in this paper employ 
the method described in Guhathakurta & Draine (1989); Sieben- 
morgen et al. (1992); Manske & Henning (1998) and Draine & 
Li (2001). For ease of reference, this section summarizes the 
method using the quantities and notation introduced in the previ¬ 
ous sections of this paper. We focus on a single grain size bin and 


where N is the number of temperature bins. Because the matrix 
values for f < i - 1 are zero these equations can be solved by 
a recursive procedure of computational order 0(N 2 ). To avoid 
numerical instabilities caused by the negative diagonal elements, 
the procedure employs a well-chosen linear combination of the 
original equations. This leads to the following recursion relations 
for the adjusted matrix elements B/, 


B ,\'~i j - A N _\j 

CN 

1 

2; 

o' 

II 

(28) 

Pfj = B f+U + A f,u 

v -s 

II 

2: 

i 

JO 

II 

o 

1 

(29) 

and for the unnormalized probability distribution A, 


Xo = 1 


(30) 

v 2}=o B ij x j 

X ' ~ Am., 

i = 1 

(31) 
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and finally for the normalized probabilities P, 

p ‘ = ^£b ■' = ».»->■ < 32 > 

^ 7=0 A 7 

While this method seems rather straightforward, specific al¬ 
gorithmic approaches differ between codes. 

One important characteristic that tends to differ between im¬ 
plementations is the grid discretizing the grain temperatures (or 
equivalently internal energy states) during the construction of 
P(T). With a fixed grid, a range of probable temperatures is de¬ 
fined, e.g. 2.7 K < T < 7 max , where 7 m;lx is chosen to exceed the 
sublimation temperature of the bulk material being considered 
and the interval is divided in N temperatures (e.g. Bianchi 2008). 
The distribution function is evaluated at that set of fixed inter¬ 
nal energy states for all grains considered to be in the stochastic 
heating regime. The advantage of this approach is that many of 
the quantities used to generate P(T) can be precomputed, reduc¬ 
ing the computational requirements of the solution method. A 
disadvantage of this approach is that calculations are done over 
the full defined temperature range, including regions where P(T) 
is negligible. Not only is this computationally inefficient, but 
it results in poor resolution of the form of P(T), especially for 
grains of intermediate size where P(T) will be relatively nar¬ 
rowly distributed with T. One alternate approach is to dynami¬ 
cally define the temperature grid (e.g. Manske & Henning 1999; 
Misselt et al. 2001). In this iterative approach, a coarse and broad 
temperature grid is defined and P(T) computed. The grid is re¬ 
fined based on P(T); temperatures with low P(T) are removed 
from the grid and P(T) is recomputed on the new, smaller, more 
densely sampled grid. The grid refinement is continued until en¬ 
ergy balance is achieved or the number of temperature points 
exceeds a predefined threshold. The advantage of this approach 
is that P(T) is properly sampled for all grain sizes. The disad¬ 
vantage is of course an increase in the computational load. This 
is amplified by the fact that the algorithm will naturally increase 
the number of temperature samples as the grain approaches the 
equilibrium regime, because P(T) is increasingly peaked as the 
grain size increases. 

A second important characteristic that tends to differ be¬ 
tween implementations is the mechanism to transition from 
stochastic to equilibrium heating regime. The simplest approach 
is to fix the grain size of the transition so that all grains with 
a < t r ans are considered to be stochastically heated and all those 
with a > a trans are considered to be in the equilibrium regime, re¬ 
gardless of the true state of the grain. Since the appropriate tran¬ 
sition point is a function of the radiation environment in addition 
to composition, this approach leads to errors in the treatment of 
the emission. However, with judicious selection of a^ans , e.g. 
a nans ~ 0.01 n m (Bianchi 2008), the results can be acceptable at 
least for non-extreme field strengths. Alternatively, the charac¬ 
teristics of P(T) can be used to terminate the stochastic heating 
treatment and transition to equilibrium heating. For example, in 
the case of a dynamic temperature grid as described above, the 
size at which the heating algorithm fails can be defined as a tran s 
(Misselt et al. 2001). A third method to determine a t rans is to 
compute the absorption and radiative timescales for each grain 
size in the considered radiation field (Draine & Li 2001). These 
timescales are a natural physical metric as stochastic heating oc¬ 
curs when the mean time between photon absorptions is long 
compared to the time the grain takes to radiatively cool. The 
ensemble of grains will then be found at a large range of tem¬ 
peratures resulting in a broad probability distribution. With this 
approach, if the absorption timescale is significantly shorter than 


the radiative timescale at a given grain size, the stochastic heat¬ 
ing regime is terminated for that and all larger sizes. 

A third characteristic that may differ between implementa¬ 
tions is the discretization of the grain size distribution, as dis¬ 
cussed in Sect. 4.1. 

5. Reference code and participating codes 

5.1. DustEM 

The DustEM code is described in Compiegne et al. (2011). 
DustEM is a stand-alone code (i.e. it is not a RT code) that 
calculates the emission and extinction of dust grains given their 
size distribution and their optical and thermal properties. It de¬ 
termines the grain temperature distribution P(T) using the for¬ 
malism of Desert et al. (1986), and it then computes the dust 
SED and associated extinction for given dust types and size 
distributions. To correctly describe the dust emission at long 
wavelengths the original algorithm has been adapted to better 
cover the low temperature region. Using an adaptive tempera¬ 
ture grid, DustEM iteratively solves the integral equation Eq. 25 
from Desert et al. (1986) in the approximation where the grain 
cooling is fully continuous. The temperature distribution calcu¬ 
lation is performed for all grain populations and sizes including 
those for which the thermal equilibrium approximation would 
apply. 

We produced the reference solutions in this paper with the 
public version of DustEM (v3.8, dated Spring 2010). To this end, 
we converted the dust properties defined in Sect. 2 and the input 
radiation fields defined in Sect. 3 into the data format expected 
by DustEM. We adjusted the values of the physical constants 
in the DustEM code, raised the maximum number of grain size 
and temperature bins to accommodate our input data, and fixed 
a minor issue in the routine that imports the grain size distribu¬ 
tion. More importantly, to obtain accurate reference solutions, 
we substantially increased the number of temperature bins and 
the number of numerical iterations (see Sect. 6.2). Other than 
this, the DustEM code was used without modifications. 

Our use of DustEM for producing reference solutions should 
not be understood to imply that it necessarily produces the phys¬ 
ically correct results. The DustEM implementation relies on the 
continuous cooling assumption just like the RT codes participat¬ 
ing in this benchmark. However, since it is not a RT code by 
itself, DustEM’s focus is solely on calculating dust properties 
and emission, and it has a neutral status in the context of this 
benchmark. 

The following sections describe each of the participating 
codes, with a focus on the specific heuristics employed for cal¬ 
culating the results presented in this paper. Table 3 offers a (very) 
concise overview of this information. 

5.2. SKIRT 

SKIRT (Baes et al. 2003, 2011; Camps & Baes 2015) is a Monte 
Carlo continuum RT code for simulating the effect of dust on 
radiation in astrophysical systems. It offers full treatment of ab¬ 
sorption and multiple anisotropic scattering by the dust, com¬ 
putes the temperature distribution of the dust and the thermal 
dust re-emission self-consistently, and supports stochastic heat¬ 
ing of small grains using an efficient library approach. The code 
handles multiple dust mixtures and arbitrary 3D geometries for 
radiation sources and dust populations, and offers a variety of 
simulated instruments for measuring the radiation field from any 
angle. It features a wide range of built-in components that can 
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Table 3. Overview of the discretization parameters and heuristics used by participating codes. Refer to Sects. 5.2 through 5.7 for a more extensive 
description. Also, Fig. 9 further supplements the information in the last column of this table. 


Code 

Grain size bins 
(Sil/Gra/PAH) 

Temperature 

bins 

Heuristic to select or determine 
temperature grid 

Heuristic to transition from 
stochastic to equilibrium regime 

SKIRT 

15/15/15 

20/625/1250 

one of 3 grids based on width of P(T) 

based on width of P(T) 

DIRTY 

121/121/28 

50-1000 

iterative range & resolution adjustment 

based on r abs /r rad 

TRADING 

20/20/8 

80 

fixed predefined grid 

fixed at a trans = 0.05 pm 

CRT 

15/15/15 

128 

one of 6 grids based on 7’ abs 

based on !P abs 

MCFOST 

63/62/24 

300 

fixed predefined grid 

based on r abs /T rad 

DART-Ray 

63/62/24 

200 

iterative range adjustment 

based on cr r of Gaussian approx. 


be configured to construct complex models without changing or 
adding source code. 

SKIRT closely follows the method presented in Sect. 4 to cal¬ 
culate dust emission. The computation time for the emissivity of 
a stochastically heated dust mixture exposed to a particular radi¬ 
ation field scales roughly with BN 2 , where B is the total number 
of size bins in the dust mixture (for all grain types combined), 
and N is the number of temperature bins used in the calculation. 

For the results shown in this paper, we use 15 size bins b for 
each grain type k, distributed logarithmically over the complete 
size range, so that B = 45. The absorption efficiencies loaded 
from the tables described in Sect. 2.4 are interpolated logarith¬ 
mically as needed to perform the integrations over grain size 
presented in Sect. 4.1 over a logarithmic grain size grid with 201 
points within each bin. 

Because of the N 2 dependency in the computation time, the 
choice of the temperature grid is rather crucial. By varying A in a 
number of experiments, it can be easily shown that, for most in¬ 
put fields, the temperature probability distribution P(T) for very 
small grains can be calculated accurately on a rather coarse grid. 
Larger grains require a finer grid because P(T') is narrower and 
has steep flanks. As discussed in Sect. 4.2, for sufficiently large 
grains P(T) approaches a delta function, so that the procedure 
described in Sect. 4.3 requires an exceedingly refined grid with 
N > 5000 to produce accurate results, which becomes compu¬ 
tationally prohibitive. Thus, in the interest of both speed and ac¬ 
curacy, we need to switch from transient to equilibrium calcula¬ 
tions for grains that can be considered to be in LTE. 

A RT simulation typically calculates the emissivity of a 
certain (fixed) dust mixture for a large number of radiation 
fields. This computation can be accelerated substantially by pre¬ 
calculating and storing the elements of the matrix A f j defined in 
Sect. 4.3, insofar as they don’t depend on the radiation field. The 
memory requirements scale with BN 2 , just like the computation 
time. More importantly, pre-calculating these values requires a 
predefined temperature grid that remains fixed for all emissivity 
calculations. However, performing all calculations on the finest 
grid would be very inefficient. 

The SKIRT implementation handles these conflicting re¬ 
quirements as follows. We predefine three separate temperature 
grids (A, B, and C) that can be used for any of the emissivity 
calculations, and we pre-calculate and store all radiation-field- 
independent values on each of these three grids for each size 
bin. The temperature range is the same for all grids; it is usually 
set to [2 K, 3000 K] but if needed the upper limit is decreased to 
the largest temperature for which enthalpy data is available in 
the dust properties. 

Grid A has only 20 bins. The widths are distributed accord¬ 
ing to a power law, providing a lot more resolution at low tem¬ 
peratures where most of the action is. The ratio of the largest 


bin (at high temperatures) over the smallest bin (at low temper¬ 
atures) is set to 500. This grid is used to find a quick estimate 
of the range in which the temperature probability distribution is 
nonzero (or rather, larger than a very small fraction of the maxi¬ 
mum probability). 

All bins in grid B are 4K wide. This medium-resolution grid 
is used to calculate the temperature probability distribution for 
dust grains with a very wide temperature range (i.e. very small 
grains and essentially all PAH grains). 

Grid C has an average bin width of 2 K, with a power-law 
ratio of 3 between the largest and smallest bins. This provides 
a fine resolution of less than 1 K at low temperatures while 
still offering decent resolution at high temperatures. This high- 
resolution grid is used to calculate the temperature probability 
distribution for dust grains with a rather narrow temperature 
range (but not so narrow that they would be considered to be 
in equilibrium). 

SKIRT implements the following heuristic to select the ap¬ 
propriate calculation for each representative dust grain: 

1. calculate the equilibrium temperature T eq for this grain; 

2. use grid A to estimate the temperature range AT = 7 max - 
T m i n in which P(T)/P mm > 1(T 20 ; 

3. if AT < 10 K or if T max < T eq then calculate the emissivity 
assuming equilibrium at temperature 7’ eq and exit; 

4. calculate P(T) using grid B (if AT > 200 K) or grid C (if 
AT < 200 K); 

5. update the temperature range AT = 7 max - T mm in which 
P(7’)/P max > 10 20 based on the new calculation; 

6. if AT ^ 10 K. oi if 7 max / c q then calculate the emissivity 
assuming equilibrium at temperature T eq and exit; 

7. calculate the emissivity using P(T) from step 4 over the 
range [T m j n , 7 nlax ] determined in step 5. 

The conditions in steps 3 and 6 are designed to avoid numer¬ 
ical instabilities when the temperature probability distribution 
approaches a delta function (AT < 10 K), and to capture situ¬ 
ations where the result is clearly inaccurate since the equilib¬ 
rium temperature lies outside the calculated temperature range 
(T max < T eq ). 

Further experiments with the SKIRT implementation show 
that for AT > 25 K the result is highly sensitive to the exact 
value of AT ; for smaller values the result converges to a sta¬ 
ble solution. For values down to AT ~ 10 K, the result is nu¬ 
merically stable in the sense that performing the calculation on 
higher-resolution grids produces essentially the same solution. 

5.3. DIRTY 

DIRTY (Gordon et al. 2001 ; Misselt et al. 2001) is a Monte Carlo 
RT code designed to study dust and its effect on radiation in arbi¬ 
trary astrophysical systems. DIRTY is a fully 3D code allowing 
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for the specification of arbitrary density distributions of both dust 
and radiation sources. It implements an adaptive mesh allowing 
for the efficient allocation of computing resources amongst re¬ 
gions in the model space depending on the physical characteris¬ 
tics of the system. Dust absorption, temperature distribution, and 
emission are handled self-consistently and multiple, anisotropic 
scattering is implemented. The dust heating implementation sup¬ 
ports both equilibrium and stochastic processes based on the lo¬ 
cal radiation field and dust properties at each grid in the model 
space. 

Like other codes presented here, DIRTY follows the ap¬ 
proach presented in Sect. 4. Internal to DIRTY, the dust grain 
size distribution is not further discretized beyond the input dis¬ 
cretization of the model; for the benchmark dust model, we com¬ 
pute the heating and emission for all sizes in the input mesh 
(28 for PAH, 121 for graphitic and silicate components; see 
Sect. 2.4). 

The heuristics employed by DIRTY in calculating the dust 
emission from each grain size of each component exposed to the 
local radiation field at a point in the model space are as follows: 

1. The equilibrium temperature, T e q , cooling timescale, r ra( j, 
and heating time scale, r a b s are computed according to sec¬ 
tion 7 of Draine & Li (2001); 

2. If the time scales computed in step 1 satisfy the inequal¬ 
ity Tabs < T ral j, the grain is considered to be in equilibrium 
with the local field; it is assigned a temperature of T eq and its 
emission is calculated following Eq. 18 with P = 5{T - T eq ). 
All grains of the same composition larger than the size for 
which the inequality is first satisfied are by default treated as 
being in equilibrium with the local field. 

3. For those grains found to be in stochastic regime in step 2, 
P(T) is computed on an initial coarse temperature grid. The 
coarse grid is defined on 50 points linearly spaced between 
[0.3,3.0] Teq. 

4. Depending on the results of step 3, the algorithm proceeds in 
one of two directions: 

(a) If P(T) is not below a specified tolerance at the endpoints 
of the temperature grid (P(7’ m j n ), P(T mix ) < lO -13 ), the 
temperature limits on the grid are expanded by 50% and 
we return to step 3 with the new, expanded temperature 
grid. 

(b) If P(T) is below the tolerance at the endpoints of the tem¬ 
perature grid, compute the energy emitted by the stochas¬ 
tically heated grain using Eq. 18. If the emitted energy is 
within 1% of the energy absorbed from the radiation field 
by the grain, consider the calculation converged and re¬ 
turn the emitted energy spectrum. Otherwise, proceed to 
step 5. 

5. The grid is now refined by a series of moves that refine the 
temperature limits and increase the number of temperature 
samples if necessary. 

(a) Remove all points on the ends of the temperature range 
which have small probabilities (P(T) < 1CL 15 ), 

(b) Recompute the probabilities on the smaller grid with the 
same number of samples. 

(c) Compute the emission and emitted energy. If the emit¬ 
ted energy matches the absorbed energy to 1%, consider 
the calculation converged and return the emitted energy 
spectrum. Otherwise; 

(d) Keeping the same temperature limits, increase the num¬ 
ber of samples by 50% and recompute the probabilities. 
Trim temperature endpoints with P(T) < 10 1 5 and re¬ 
turn to step 5b. If this step would result in the number 


of bins exceeding a pre-defined maximum (1000), record 
the failure of the stochastic heating algorithm and pro¬ 
ceed to the next grid in the model space. 

In practice, the failure described in step 5d occurs in model 
bins for which the local field has not been well defined, generally 
due to very few photons interacting in that cell, either through 
a poor definition of the adaptive mesh or insufficient photons 
being run in the Monte Carlo simulation of the radiation. These 
cells are generally unimportant in the overall energy budget of 
the model and can be masked in post processing. Such cells are 
rare in most model runs; their number and distribution in the 
model space can be used as a metric of the overall quality of the 
simulation. 

The approach in step 2 results in all PAH sizes being treated 
in the stochastic regime in the radiation fields explored in this 
benchmark. Silicate and graphitic grains generally achieve equi¬ 
librium between sizes of 0.006 and 0.020 /urn. 

5.4. TRADING 

TRADING (Bianchi 2008) is a 3D Monte Carlo dust contin¬ 
uum RT code with characteristics similar to those of SKIRT and 
DIRTY. Originally designed to study the effects of clumping in 
the disks of spiral galaxies, it uses an binary-tree adaptive grid 
(octree) for the dust distribution. 

TRADING computes stochastic heating following the 
method described in Sect. 4. A single temperature grid is used 
to precompute the field-independent terms of the matrix ele¬ 
ments (Eq. 25 and most factors in Eq. 23). Since the RT models 
of clumpy galactic disks in Bianchi (2008) need a few million 
dust cells, and each cell requires the calculation of dust heat¬ 
ing for the full grain size distribution, Bianchi (2008) uses a 
limited temperature grid of 80 logarithmically spaced bins be¬ 
tween 2.7 and 2000 K. As mentioned in Sect. 4.3, faster thermal 
equilibrium calculations were performed for grains larger than 
fltrans = 0.01 /mi. The original setup resulted in SEDs that were 
estimated to be within 10% of full solutions for wavelengths up 
to 1000 fjm. 

For this benchmark we left the number of temperature bins 
at 80, but we extended the temperature range to 3000 K. While 
the previous choice of atmns produces accurate results for the 
typical interstellar radiation fields encountered in spiral galax¬ 
ies (with U > 0.1, see Aniano et al. (2012); Hunt et al. (2014)), 
we adopted here a lrans = 0.05 fiva in order to improve the solution 
for fields with lower intensities. 

The dust grain size distribution was discretized using a grid 
with 20 size bins for graphite, 20 for silicates, and 8 for PAHs, 
logarithmically spaced over their size range. The choice allows 
to have similar bin widths for all materials. Optical properties 
were derived interpolating logarithmically over those of the full 
size table. The grain size grid used for this paper is similar to 
that adopted in Bianchi (2008) though for a different dust model. 

For a t rans = 0.01 //m, only about half of the size bins pass 
through the stochastic heating calculation. For a trans = 0.05 //m, 
this goes up to 75%: all of the PAH bins and 14 of the 20 graphite 
and silicate bins. If a grain with a < a tran s attains the condition 
of thermal equilibrium, the adopted temperature grid might fail 
in computing the resulting narrow probability function. Thus, in 
addition to the grain size cut-off, we chose to assume thermal 
equilibrium whenever the range of the computed temperature 
distribution defined by P(T) > 10 1 5 is smaller than the equi¬ 
librium temperature. 
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5.5. CRT 

CRT (Juvela & Padoan 2003; Juvela 2005; Lunttila & Juvela 
2012) is a Monte Carlo dust continuum RT program. It solves 
the RT equation self-consistently with a full treatment of scat¬ 
tering, absorption and emission of radiation in ID, 2D, and 3D 
geometries. The program allows using an external component, 
for instance DustEM, for the calculation of the dust emission 
spectrum for a given input radiation field. This feature has been 
used, e.g. in Ysard et al. (2011), to study the microwave emis¬ 
sion from spinning dust grains. However, CRT itself also has 
optimized routines for fast dust emission calculations, including 
the treatment of stochastically heated grains. Although CRT can 
calculate dust emission with fully discrete cooling, the results 
presented in this paper use the continuous cooling approxima¬ 
tion, and the following discussion is focused on the algorithms 
for continuous cooling computations. 

The basic algorithms employed by CRT follow the outline 
presented in Sect. 4. To allow the use of pre-computation to 
speed up the construction of the transition matrix A , the tem¬ 
perature grid is not defined fully dynamically according to the 
input field. Instead, each dust type and size uses one of several 
predefined grids, for which the pre-computations are done at the 
beginning of calculation. The predefined grids are linear in T and 
their upper and lower limits are chosen to allow a good represen¬ 
tation of the grain temperature distribution in the types of radi¬ 
ation fields that are found in the model. In particular, the grids 
should be built for reference fields that span approximately the 
range of radiation energy densities expected to be found in the 
model. If the hardness of the radiation field varies significantly, it 
may also be useful to include grids for different spectra with the 
same energy density. For the calculations presented in this paper 
we use six temperature grids that were built for scaled Mathis 
ISRF with U = ICC 4 , KT 1 ,1,10,10 3 , and 10 6 . For large grains 
and strong radiation fields, the predefined grid has a special en¬ 
try that triggers equilibrium calculation, otherwise full stochastic 
calculations is used. The number of temperature bins for stochas¬ 
tic calculation can be set by the user; in this paper we use 128 
bins. 

To select the temperature grid that is used for calculating 
the emission for a grain size and type in a given input field, we 
calculate the absorption time scale r a b S and the mean energy of 
the absorbed photon (£abs)- We choose the predefined grid that 
has been built for the reference field most like the current input 
field. In the calculations presented in this paper the selected grid 
is simply the one for whose reference field the mean absorbed 
power per grain P^bs = (Labs) / T abs is closest to the values calcu¬ 
lated for the input field. Although the selected temperature grid 
is not necessarily optimal, it is good enough to allow accurate 
results using a modest number of energy bins. 

The computation of emission from stochastically heated 
grains in CRT differs slightly from the description given in 
Sect. 4.3. Instead of using Eq. 23 for calculating the upwards 
transition rates, we apply Eqs. 15-25 and 28 from Draine & 
Li (2001), which include corrections for the finite size of an 
enthalpy bin. Similarly, instead of Eq. 25, we use Eq. 41 from 
Draine & Li (2001) for calculating the cooling part of the tran¬ 
sition matrix. Including these finite bin size corrections allows 
using a lower resolution grid, which substantially benefits the 
computation time. 

The cooling matrix elements (/ < i) are independent of the 
radiation field and they are pre-computed. The heating elements 
(/ > i ) depend on the radiation field and must be calculated sep¬ 
arately for each input field. Moreover, when using the equations 


from Draine & Li (2001), instead of using the radiation field 
strength at a single wavelength as in Eq. 23, we must integrate 
numerically over the wavelength grid. We use precomputed in¬ 
tegration weights Wfj corresponding to each grid point of the 
wavelength grid. The integral in Eq. 15 in Draine & Li (2001) 
can then be evaluated as Afj - w f,i(k)J(Ak ). If the number 

of points in the wavelength grid is large, calculating the full sum 
is slow. However, for given energy bins i and /, radiation within 
only a narrow range of wavelengths can induce transitions i —> f. 
Therefore, Wf j(k) is non-zero only for a few k and only non-zero 
integration weights are stored and used in the summation. 

Discretization of the particle size distribution can be defined 
by the user. In this paper we employ the same discretization 
as SKIRT, i.e., 15 logarithmic bins for each of the three dust 
components. Dust properties for each size bin are calculated ac¬ 
cording to Eqs. 15-17 using numerical integration with 256 grid 
points. 

5.6. MCFOST 

MCFOST (Pinte et al. 2006, 2009) is a 3D continuum and line 
RT code. It relies on the Monte Carlo method to compute the lo¬ 
cal specific intensities and related quantities (e.g. temperature, 
molecular levels) and computes observables via a ray-tracing 
method. The emerging fluxes are calculated by formally inte¬ 
grating the RT method along rays using the specific intensities 
and source functions computed during the Monte Carlo run. 

MCFOST computes the stochastic heating of small dust 
grains following the method presented in section 4, with a few 
refinements to ensure numerical stability and speed. 

We first compute the time between two successive absorp¬ 
tions of a photon and compare it to the cooling time of the grain, 
following the method described in Draine & Li (2001). For dust 
grains where the time between two absorptions is smaller than 
the cooling time, only the equilibrium temperature is calculated. 

For those grains which have a shorter cooling time, we com¬ 
pute the full temperature probability distribution using a fixed 
temperature grid with 300 points logarithmically distributed be¬ 
tween 1 and 3000 K. The cooling terms of the transition matrix 
(Eq. 25), which are independent of the radiation field, are pre¬ 
computed, as well as most of the heating terms factors (Eq. 23). 
We estimate the specific intensity J(Afi) for each term in Eq. 23 
by interpolating the radiation field computed by the Monte Carlo 
run. As the interpolation coefficients are identical for every cell 
in the model, they are also precomputed to speed-up the calcula¬ 
tions. 

For dust in radiative equilibrium, MCFOST solves the prob¬ 
lem of self-consistent dust heating and re-emission using the im¬ 
mediate re-emission concept of Bjorkman & Wood (2001). This 
methods eliminates the need for iteration and ensures a perfect 
conservation of energy. For non-equilibrium dust grains how¬ 
ever, this procedure is prohibitive as it requires a temperature 
calculation at each absorption/re-emission. Instead, we use the 
classical iterative scheme where we store the energy absorbed by 
the dust, compute the temperature probability distribution in all 
cells once all the packets have been propagated and re-emit the 
absorbed energy via new packets according to the new temper¬ 
ature probability. The procedure is iterated until a desired con¬ 
vergence on the temperature or energy is reached. Because only 
the fraction of radiation that is absorbed by the non-equilibrium 
grains needs to be re-emitted in an iterative way, convergence 
is usually reached after only a few iterations. In practice, when 
a packet is absorbed inside a cell, its energy is split: the frac- 
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tion absorbed by dust grains in equilibrium is immediately re¬ 
emitted, while the fraction absorbed by non-equilibrium grains 
is stored to be re-emitted during the next iteration. For the bench¬ 
mark presented in this paper, the input radiation field is fixed so 
no iteration is required. 

In the presented calculations, the grain size distributions 
were discretized using 63 logarithmically spaced grain sizes for 
silicates, 62 for graphite and 24 for PAHs. 

5.7. DART-Ray 

DART-Ray is a ray-tracing 3D dust RT code that implements the 
RT algorithm described in Natale et al. (2014). It can be used to 
derive radiation field energy density distributions and outgoing 
radiation surface brightness maps for arbitrary 3D distributions 
of dust mass and stellar emission. It includes treatment of both 
absorption and anisotropic scattering. For the dust emission cal¬ 
culations, DART-Ray uses the prescription initially incorporated 
in the 2D RT model of Popescu et al. (2000), and later updated in 
Popescu et al. (2011). However, unlike the 2D models where the 
stochastically heated dust emission could be explicitly computed 
for each individual position, the calculations for the 3D models, 
containing of the order of 10 6 cells, can be accelerated by using 
an adaptive SED library approach (see Natale et al. 2014b, in 
press). 

For a given radiation field intensity spectrum, found for a 
particular model cell, the stochastically heated dust emission is 
derived following Voit (1991). The method used to determine the 
probability distribution P(T) combines the numerical integration 
of Guhathakurta & Draine (1989) with a step-wise analytical so¬ 
lution. The algorithm provides accurate and swift results on a 
relatively coarse grid. This is particularly useful for larger grain 
sizes where the probability distribution P(T) converges to a nar¬ 
row distribution around the equilibrium temperature. In case of 
the first order integration of Guhathakurta & Draine (1989), the 
width of the energy bins needs to be considerably smaller than 
the mean deposited energy to preserve energy balance. An in¬ 
creasing number of energy bins would be required to avoid en¬ 
ergy losses in the heating process, which would make the calcu¬ 
lation inefficient (see Fischera 2000, and Sect. 4.3). 

The absorption of CMB photons is assumed to provide a con¬ 
tinuous heating source. It is taken into consideration by subtract¬ 
ing the heating rate related to the CMB from the cooling rate 
(Eq. 25), which limits the temperatures to values not lower than 
the CMB temperature. As further discussed in Sects. 6.3 and 6.6, 
the cooling assumed in the dust model is only valid as long as 
the emitted photon energy of the modified black body is small 
relative to the enthalpy of the grain. This assumption does no 
longer apply at very low temperatures of small dust grains or 
PAH molecules (Fischera 2000). While negligible for most cases 
we find that the CMB becomes a considerable heating source for 
the very low radiation field strength. For U = 10 4 the heating of 
silicate grains by CMB photons is for all sizes larger than 10%. 
For graphites the contribution is considerably lower with only a 
few percent and still negligible for PAH molecules. 

The temperature distribution P{T) for each grain size of a 
given composition is obtained consecutively by starting with the 
largest grain size and utilizing the basic characteristic that the 
distributions broaden with decreasing grain size. The probabil¬ 
ity distribution is only considered for all grains below a certain 
size where the stochastic heating leads to a considerable distri¬ 
bution of the dust temperatures. Above this critical grain size the 
grains are assumed to radiate at the equilibrium temperature T eq . 
To estimate the transition from equilibrium to non equilibrium. 


we apply the Gaussian approximation in the limit of large grains 
as derived by Voit (1991). We consider the grains to be stochas¬ 
tically heated if Ictt/T eq > 0.05. For 0.05 < 2cr T /T eq < 0.1 we 
apply the Gaussian approximation; for 2cr T /T eq >0.1 the full 
P{T) distributions are derived. 

The temperature distributions of stochastically heated dust 
grains are derived for dynamically determined temperature inter¬ 
vals [Tmin, 7 milx ]. For the results presented in this paper, we sub¬ 
divide the temperature interval into 200 temperature bins equally 
spaced in log T. The interval boundaries are determined itera¬ 
tively by increasing T m i n or decreasing 7 max by 30% until the 
probabilities P(T) for the lowest and highest temperature bin are 
lower than 10 20 , ensuring that the emission at higher and lower 
temperatures can be neglected. To accelerate the iterative pro¬ 
cess, we use for each dust grain as the initial estimate of the tem¬ 
perature interval the derived interval of the previous larger dust 
grain. For the first grain for which the temperature distribution is 
derived, we use as initial estimate of the temperature interval a 
width based on the Gaussian approximation in the limit of large 
dust grains. 

For the results presented in this paper, we derive the dust 
emission for each grain species at each grain size of the tabu¬ 
lated size distribution described in Sect. 2.4. To calculate the to¬ 
tal emission, we then integrate over the size distribution and sum 
the contributions from different grain species. By comparing the 
total dust emission with the total absorbed energy we ascertain 
that the energy balance for every grain species is fulfilled with 
an accuracy better than a few percent. The largest discrepancies 
(3-4%) are found for the most extremely scaled Mathis fields 
considered in the benchmark (U = 10 4 and U = 10 6 ). 

6. Results and discussion 

6.1. Data files 

The data files representing the benchmark results can be down¬ 
loaded from the web site indicated in footnote 2. For each par¬ 
ticipating code, and for each input radiation field, the calculated 
solution is stored in a separate text file with columns specifying 
the wavelength A (in /vrn); the mean intensity of the input field 
(in W m 3 sr 1 ); and the silicate, graphite and PAH emissivities 
As,i (in W sr -1 H 1 ), in that order. The file naming scheme and 
the precise file format are described on the web site. 

6.2. Reference solutions 

As mentioned in Sect. 5.1, below we compare the results from 
the codes participating in this benchmark against reference so¬ 
lutions generated with DustEM. Figure 1 shows these solutions 
for a selection of the input fields defined in Sect. 3. 

To ensure proper accuracy, we increased the number of tem¬ 
perature bins and the number of iterations in the DustEM integral 
equation solver until the calculated emission converged to a sta¬ 
ble solution; see Fig. 2. Specifically, the number of temperature 
bins was raised from 200 to 3500, and the number of iterations 
from 80 to 250. These changes dramatically increase the com¬ 
putation time, however this is acceptable for the calculation of a 
reference solution. 

As an extra sanity check, we verified that the emissivities 
calculated by DustEM (using the same number of tempera¬ 
ture bins and iterations as for the reference solutions) indeed 
converge to the corresponding equilibrium emissivities. Fig. 3 
shows this comparison for 0.05 /rm grains exposed to radiation 
fields ranging from extremely weak (left) to strong (right). For 
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a strong field, where we expect the grains to be in equilibrium, 
the DustEM solutions indeed match the equilibrium emissivi- 
ties. For a weak field, the solutions differ since the grains are 
no longer completely in equilibrium. This shows the importance 
of performing the full stochastic calculation in the presence of 
extremely weak fields, even for grain sizes up to 0.05 gm. Codes 
transitioning from stochastic to equilibrium calculation at a fixed 
grain size should thus set a sufficiently high value of a ttans (see 
Sects. 4.3 and 5.4). 

6.3. Benchmark solutions 

Fig. 4 compares the total emissivities calculated by each of the 
codes participating in this benchmark to the corresponding ref¬ 
erence solutions for a selection of the input fields defined in 
Sect. 3. Subsequent figures zoom in on the emissivities for each 
dust component separately, i.e. silicate (Fig. 5), graphite (Fig. 6), 
and PAH (Fig. 7) grains. 

In these figures we limit the displayed wavelength range to 
1 pm < A < 1000pm. Outside of this range, other sources or 
processes usually dominate the radiation emanating from astro- 
physical objects, so it is not relevant to evaluate the results of 
a dust emissivity calculation in that spectral range. Moreover, 
some of the assumptions underlying the computations are no 
longer valid, rendering the results physically meaningless. 

First considering the shorter wavelength range, a black body 
with peak emission at A - 1 pm has a temperature of T ~ 
2900 K. The sublimation temperature of a dust grain is estimated 
at 1200K for silicates and at 2100K for graphites (Kobayashi 
et al. 2009). Evaporation rates rise roughly exponentially with 
increasing grain temperature (Guhathakurta & Draine 1989; 
Kobayashi et al. 2009), i.e. with decreasing emission wave¬ 
length. It is clear that a relevant fraction of the dust that would 
emit at /l < 1 pm is destroyed by evaporation, so that the grain 
size distribution in our dust model is no longer valid under these 
conditions. Consequently, the calculation would substantially 
overestimate the resulting dust emission. 

We now consider the longer wavelength range. For a suffi¬ 
ciently strong input field, say U > 1, we can expect the calcu¬ 
lated emissivity results to be correct because most of the grains 
emitting at wavelengths longer than 1000/mi are in LTE. How¬ 
ever, the emissivity peaks at much shorter wavelengths, so that 
the level at 1000 pm is already several orders of magnitude be¬ 
low the peak level (see all panels in Fig. 1 except for the first 
two). The situation is different for extremely weak input fields 
approaching the level of the cosmic microwave background 
(CMB). A black body with peak emission at A = 1000 pm has 
a temperature of T ~ 2.9 K, just above the temperature of the 
CMB. For small dust grains at such low energies, the continuous 
cooling assumption no longer holds, 4 meaning that the method 
presented in Sect. 4.3 does not necessarily yield correct results. 
In conclusion, the dust emission at wavelengths A > 1000 pm 
is either calculated assuming equilibrium conditions (which ren¬ 
ders comparison uninteresting) or calculated improperly (which 
renders comparison meaningless). 

6.4. Evaluation of benchmark results 

In the wavelength range 3 yum < A < 1000 pm all participat¬ 
ing codes reproduce the total dust emissivity within 20% of the 

4 With the properties of our dust model, the internal energy of a 
0.007 pm graphite grain at 2.9 K is insufficient to emit a 1000 pm pho¬ 
ton. 


reference solution for all input fields used in this benchmark 
(see Fig. 4). Excluding the weakest (U < 10 4 ) and the soft¬ 
est (T < 3000 K) fields, the correspondence in the same wave¬ 
length range is within 10%. The larger relative deviations at 
wavelengths shorter than 3 pm are caused in part by the much 
lower absolute emissivity values in that range (two to three or¬ 
ders of magnitude below peak values; see Fig. 1 and the vertical 
dashed lines in Fig. 4). 

The emissivities calculated for the silicate and graphite com¬ 
ponents (Figs. 5 and 6) show a similar pattern, although the de¬ 
viations in the individual components are sometimes slightly 
larger. The emissivities calculated for the PAH component 
(Fig. 7) show the largest deviations. If we restrict the analysis 
to the wavelength range in which the emissivity of the reference 
solution is within three orders of magnitude of its peak value (in¬ 
dicated by the vertical dashed lines), the correspondence is still 
within 40% for all codes and for all input fields, and often a lot 
better. 

The larger discrepancies between the various codes for the 
PAH component can be traced to the fact that PAH molecules 
are, generally, substantially smaller than silicate and graphite 
grains; see the a"' ax values in Table 2. As noted in Sect. 4.3 and 
further discussed in Sect. 6.5, smaller grains are more likely to 
remain in the stochastic regime, requiring complex calculations 
that are more sensitive to differences in discretization (choice 
of grids) and concrete implementation (even if the same overall 
method is employed), as compared to equilibrium calculations. 

Specifically, the PAH emissivities calculated by CRT (ma¬ 
genta curves in Fig. 7) deviate from the other codes because, as 
described in Sect. 5.5, CRT implements the Draine & Li (2001) 
equations rather than Eqs.23 and 25 for calculating the heat¬ 
ing and cooling transition rates. This method does not allow the 
emission of photons with an energy higher than the enthalpy con¬ 
tent of the emitting grain (see Eq. 56 in Draine & Li 2001). The 
upwards jumps in the emission spectrum appear when, going 
towards longer wavelengths (lower photon energy), a new en¬ 
thalpy bin enters the emission calculation (see also Figs. 14 and 
15 in Draine & Li 2001). These discontinuities appear only in 
the wavelength range where emission is largely from grains with 
very low enthalpy, for which the continuous cooling approxima¬ 
tion is not valid, as described in Sect. 6.3. 

Fig. 8 shows the contributions to the total emissivity calcu¬ 
lated by one of our codes (DIRTY) in the stochastic and equilib¬ 
rium regimes for each of the three grain types, and for a number 
of input fields. '' While the silicate and graphite components show 
a significant equilibrium contribution for all fields, the PAHs re¬ 
main in the stochastic regime for all but the strongest fields. 

With respect to our conclusion at the end of Sect. 6.3, it is 
worth noting in Fig. 8 that, at A = 1000 pm. the equilibrium 
contributions of the silicate and graphite components dominate 
the total stochastic contribution for all input fields. Only for the 
weakest fields (U ~ 10 4 ) does the stochastic contribution at that 
wavelength become a noticeable fraction of the total emission. 

Considering the contributions of the different grain types to 
the total spectrum (Fig. 1), it turns out that, for our dust model, 
the graphite emission dominates over most of the wavelength 
range for most input fields. The PAHs dominate in a very small 
region (their peak) and the silicates only at the longest wave¬ 
lengths. Since the agreement between the different codes is sig¬ 
nificantly better for graphite and silicates than for the PAHs, 


5 The precise form of these respective contributions varies between 
codes because of differences in transitioning from one regime to the 
other, but the general trend is similar. 
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our choice of dust mixture (unintentionally) benefits the global 
agreement between the different codes. 

6.5. Transition to equilibrium 

As introduced in Sect. 4.3, and further elaborated upon in the 
code descriptions in Sect. 5, each code handles the transition 
from the stochastic to the equilibrium calculation regime in its 
own way. Figure 9 shows the grain size a trans for which the par¬ 
ticipating codes transition from the stochastic to the equilibrium 
calculation regime, for each grain type, and for each of the input 
fields defined in Sect. 3. While the details differ between codes, 
as a general trend, small grains (e.g. a < 0.01 pm) are consid¬ 
ered to be in equilibrium only when exposed to the strongest 
fields (U > 10 2 ). 

Because most of the PAHs remain in the stochastic regime 
(Fig. 8 and right panel of Fig. 9), the transition differences are 
most easily seen in the results for the silicate and graphite grains 
(Figs. 5 and 6). As a general trend, the results converge for 
longer wavelengths, because the equilibrium emission dominat¬ 
ing in that range is calculated in the same straightforward man¬ 
ner across all codes. The larger discrepancies are found at shorter 
wavelengths, where the stochastic regime dominates. Depending 
on the input field, the transition point shifts on the wavelength 
scale. Interestingly, for some fields, the discrepancies show ex¬ 
tra „wiggles” near the transition points, most likely caused by 
differences in handling the transition. This is particularly evi¬ 
dent in the silicate emissivities for the scaled Mathis fields with 
strengths U = 1CF 2 to 10 2 , (see the top half of Fig. 5). 

The total emissivities are influenced by these transition dif¬ 
ferences mostly in the wavelength range just short-ward of the 
large submm emission peak, which is dominated by LTE emis¬ 
sion from silicate and graphite grains (see Fig. 4). The position 
of this peak is determined by the temperatures of the dust grains, 
and thus depends on the input field. For the strongest field in our 
benchmark (U = 10 6 ), the broad peak even overlaps the PAH 
features in the 3 to 30 pm wavelength range (see the third panel 
in the leftmost column of Fig. 4). 

6.6. Weak fields 

As discussed in Sect. 4.3, the transition size a lrdm above which a 
grain can be considered to be in LTE depends on the radiation 
field to which the grain is exposed. For a given grain composi¬ 
tion, a trans tends to be higher for weaker fields because photon 
interactions are less frequent, which keeps larger grains in the 
stochastic regime. Consequently, the dust emission calculations 
are more complex for weaker fields, and especially for small 
grains in weaker fields. In addition to the computation time, this 
complexity affects the accuracy of the results, which explains 
why the largest discrepancies between the various solutions oc¬ 
cur for the weakest field in the benchmark (see the top left panels 
in Figs. 4 through 7). 

In fact, the weakest field in our benchmark (U - 10 4 j may 
be unrealistically weak, since its peak intensity is below the peak 
of the CMB - albeit in a different wavelength regime (UV ver¬ 
sus mm wavelengths). To evaluate the effect of neglecting the 
CMB, we added the CMB to the 10~ 4 x 7 Mat field, and used 
DustEM to recalculate the emissivity of our dust model exposed 
to this new input field. For wavelengths A < 100/rm the results 
are essentially identical to those shown in the top left panel of 
Fig. 1. The submm peak is a notch higher and slightly shifted 
to longer wavelengths, causing an emissivity increase of about 


35% at A = 1000yum. While this effect may not be negligible, it 
does not invalidate the benchmark test. 

However, as argued in Sect. 6.3, our computations may no 
longer be physically founded for these weak fields, especially 
for small grains with internal energies comparable to those of the 
CMB photons. In conclusion, the 10 4 x 7 Mat input field is prop¬ 
erly benchmarking the various codes, but the calculations may 
be collectively incorrect because the continuous cooling approx¬ 
imation is inappropriate in this regime. 

6.7. Temperature discretization 

The participating codes implement various ways to discretize the 
grain temperature (or equivalently, the grain enthalpy), as de¬ 
scribed in Sect. 5. The different schemes are mostly driven by the 
aim to increase performance while preserving accuracy. Here we 
discuss the impact of respectively the minimum and maximum 
temperature values allowed in the grid. 

DustEM, used to calculate our reference solutions, does not 
impose a lower temperature limit other than the zero point. In¬ 
deed, under the continuous cooling assumption, a small dust 
grain does not have to be in equilibrium with the CMB, and thus 
there seems to be no reason why the grain should not have, at 
a given moment in time, a temperature below 2.73 K. In other 
words, we need to calculate the temperature probability distribu¬ 
tion as usual. As argued in Sects. 6.3 and 6.6, this line of reason¬ 
ing breaks down for small grains at very low energies, since the 
continuous cooling assumption no longer holds. 

This is why most codes participating in this benchmark do 
impose a lower temperature limit of 2.73 K, or even, rather ar¬ 
bitrarily, 1, 2 or 3K. A limit of 2.73 K, for example, causes a 
bump in the PAH emission peaking at A « 1060pm (the peak 
CMB wavelength) because all the probabilities for lower tem¬ 
peratures are bunched together in the 2.73 K temperature bin. 
This effect is to some extent responsible for the discrepancies 
between the codes and the reference solution seen in the top left 
panel of Fig. 7, in the wavelength range to the right of the dashed 
line, where the absolute value of the emissivity has become small 
anyway. The effect is negligible for all but the weakest fields. 

At the other end of the scale, all codes in this benchmark 
use the complete temperature range for which the dust proper¬ 
ties are defined, i.e. up to 2500 K for our dust model. As de¬ 
scribed in Sect. 6.3, this is well above the sublimation tempera¬ 
ture of the dust material, although a fraction of the grains may 
survive at these temperatures for some time. Because the method 
used for this benchmark ignores dust grain destruction, we ex¬ 
pect it to overestimate the emissivity for shorter wavelengths. To 
evaluate this effect, we reran the benchmark calculations with 
one of the codes (SKIRT) using a maximum grid temperature of 
2250 K instead of 2500 K. As expected, the emissivities in the 
range ,1 > 3 pm are essentially unaffected by this change. At 
wavelengths shorter than 3 pm, the total emissivity for the hard¬ 
est black body input fields (T > 15000 K) decreases noticeably 
(by about 30% at 1 pm), while there is no perceptible change 
down to 1 pm for the softer fields or for the scaled Mathis fields. 
The silicate and graphite components behave similarly. Interest¬ 
ingly, the emissivity of the PAH component shows a noticeable 
decrease for all but the softest black body fields (T < 6000 K), 
and for all scaled Mathis fields. This is can be understood by re¬ 
calling that the PAH particles are, on average, a lot smaller than 
those in the other components, so that they are more easily pro¬ 
pelled to higher temperatures. 
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6.8. Wavelength discretization 

The method described in Sect. 4 for calculating dust emission 
involves wavelength discretization in several distinct areas. The 
optical dust properties are tabulated on some predefined wave¬ 
length grid. We also need to configure the wavelength range and 
sampling resolution for the input fields and for the output emis- 
sivity. And finally, the calculation of the cooling coefficients de¬ 
fined in Eq. 25 requires an integration of the optical properties 
over wavelength, which implies a grid as well. As long as proper 
interpolation procedures are in place, these wavelength grids do 
not need to be identical. 

To keep matters simple for the benchmarks presented in this 
paper. Sect. 3.3 specifies the same wavelength grid for the optical 
properties and for the calculated emissivities, i.e. a logarithmic 
grid with 1201 points in the range 0.001 pm < A < 10000 pm. 
In this section we discuss the impact of the resolution and of the 
lower and upper limits of the wavelength grid on the emissiv- 
ity calculations. We reran the benchmark calculations with one 
of the codes (SKIRT) using some wavelength grid variations as 
reported below. SKIRT employs a single (configurable) wave¬ 
length grid for all aspects of the calculation. The optical dust 
properties are interpolated to this grid, and the grid is subse¬ 
quently used for the input/output fields and for the integration 
to obtain the cooling coefficients. 

In practice, the lower wavelength limit is determined by the 
need to properly capture the input field in the calculation of 
the heating coefficients (Eq. 23) and the equilibrium tempera¬ 
ture (Eq. 22). For a scaled Mathis field, the lower limit can be 
increased to 0.09 pm (see Eq. 11). For the hardest black body 
fields (T > 12000 K), the limit should be lower. With a lower 
limit of 0.01 pm there is no noticeable difference in the calcu¬ 
lated emissivities even for the hardest field in this benchmark 
(T = 18000 K) 

The upper wavelength limit is mostly determined by the need 
to calculate the emissivity up to mm wavelengths. The upper 
limit also affects the calculation of the heating coefficients and 
the equilibrium temperature, but this effect is smaller because 
the absorption coefficients for the grain material are much lower 
at longer wavelengths. With an upper limit of 2000 pm instead 
of 10000 pm there is no noticeable difference in the calculated 
emissivities for any of the input fields. With an upper limit of 
1000 pm the submm emissivity peak is overestimated by 20% 
for the weakest scaled Mathis field (U = ICE 4 ). The emissivity 
peak for the U = 10 4 field shows a similar but much smaller 
effect. For all other input fields there is no noticeable difference. 

Finally, we reran the benchmark calculations for logarithmi¬ 
cally distributed wavelength grids with successively lower reso¬ 
lution, always using a range of 0.01 pm < A < 2000 yum. First of 
all, lowering the wavelength resolution affects the shape of the 
sharp PAH dominated emissivity peaks in the 3 to 30yum wave¬ 
length range; if there is no wavelength point at the center of a 
peak, the peak can not be resolved. However, this does not af¬ 
fect the accuracy of the emissivity at other wavelength points 
(unless the resolution becomes too low, as described in what fol¬ 
lows). Other than this peak resolution effect, using 601 wave¬ 
length points instead of 1201 does not noticeably influence the 
results for any of the input fields. Fowering the resolution to 301 
wavelength points causes minor deviations in the calculated PAH 
emissivities for wavelengths A < 3 pm. which however do not 
noticeably affect the total emissivities down to A > 1 pm. With 
only 151 wavelength points the deviation in the total emissivi¬ 
ties is still limited (a few percent at 1 yum), but the PAH features 
are now clearly under-resolved and rather smoothed out. This 


can be improved by concentrating more grid points in the wave¬ 
length range of the PAH features. For example, a specialty grid 
with a total of 151 points, 61 of which are concentrated in the 
3 to 30 pm wavelength range, seems to provide an acceptable 
compromise. 

6.9. Grain size discretization 

The total calculation time for the emissivity of a dust popula¬ 
tion is roughly proportional to the number of grain size bins 
used to represent the population (see Sect. 4.1). Therefore, some 
of the participating codes recompute the optical grain proper¬ 
ties on an internally defined grain size grid rather than using the 
size bins tabulated in the dust model data. We used two of the 
codes (SKIRT and CRT) to investigate the effect of the number 
of grain size bins on the calculated emission. As mentioned in 
Sects. 5.2 and 5.5, these codes were configured with 15 size bins 
per grain type to calculate the benchmark results presented in 
Figs. 4 through 7. 

The predominant effect of changing the number of grain size 
bins appears for wavelengths A < 10/rm and grows larger for 
shorter wavelengths. This is to be expected; with a coarse grid 
the effective size of grains in the smallest bin is relatively large, 
and it is difficult to heat the grains to the high temperatures 
that are needed for emission at shorter wavelengths. Specifically, 
when the number of size bins per grain type is reduced from 15 
to 10, the calculated silicate and graphite emissivities show sub¬ 
stantial deviations from the reference solutions. The PAH emis¬ 
sivities are virtually unaffected, which is easily understood be¬ 
cause their size range is much smaller. When the number of size 
bins per grain type is increased to 30, the calculated silicate and 
graphite emissivities in the wavelength range 1 pm < A < 10 pm 
approach the reference results. Again the PAH emissivities are 
much less affected because the dust model data has only 28 PAH 
size bins anyway. 

A secondary effect occurs for longer wavelengths because 
the re-binning influences the heuristic for transitioning between 
the stochastic and equilibrium calculation regimes (Sect. 6.5). 
This effect seems to be somewhat random in nature, causing 
deviations that remain within the accuracy limits described in 
Sect. 6.4. 


6.10. Calculation time 

A typical 3D RT simulation calculates the dust emission spec¬ 
tra for a large number of dust cells. In cases where dust self¬ 
absorption is a relevant factor, this calculation is repeated for 
each iteration of the loop that self-consistently determines dust 
heating and re-emission (see e.g. Sect. 5.6). The time spent on 
calculating dust emission might thus become a significant or 
even dominant fraction of the total RT simulation time. The aim 
to reduce the dust emission calculation time has guided many of 
the choices in the implementations of the RT codes participat¬ 
ing in this benchmark. Most fundamentally, all codes adopt the 
continuous cooling approximation. In addition, all codes select 
specific discretization schemes, most codes employ heuristics to 
transition between stochastic and equilibrium regimes, and some 
pre-compute field-independent data. Often these choices affect 
not only the calculation time, but also the accuracy of the results. 
In principle at least, the results can be made to match perfectly 
by increasing grid sizes and removing the heuristics, at the ex¬ 
pense of calculation time. 
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It thus seems appropriate to consider calculation time when 
evaluating benchmark results. Unfortunately, it is not meaning¬ 
ful to compare the dust emission calculation times between the 
codes outside of the context of a RT simulation. For example, 
moving the relevant data for each dust cell from memory into 
the processor cache and back may represent a significant por¬ 
tion of the total calculation time, depending on the memory lay¬ 
out chosen by the RT code and depending on the architecture of 
the computer system. Consequently, a performance comparison 
would be more appropriately conducted as part of a RT bench¬ 
mark. 

Just to provide an order of magnitude, with the prescriptions 
provided in this paper, a code can calculate a few hundred dust 
emission spectra per second on a modern desktop computer. This 
means that the calculation for 5 million dust cells can be com¬ 
pleted in a matter of hours rather than days. 

7. Conclusions 

We defined an appropriate problem for benchmarking dust emis- 
sivity calculations in the context of RT simulations, specifically 
including the emission from stochastically heated dust grains 
(SHGs). The problem definition includes the optical and calori¬ 
metric material properties, and the grain size distributions, for 
a typical astronomical dust mixture with silicate, graphite and 
PAH components; a series of analytically defined radiation fields 
to which the dust population is to be exposed; and instructions 
for the desired output. 

We summarized a popular method for calculating the emis¬ 
sion from SHGs, with the intention to provide a self-contained 
guide for implementors of such functionality. The method is fre¬ 
quently used in RT codes because of its good performance and 
relative ease of implementation, although it assumes continuous 
cooling of the dust grains, which may be inaccurate in extreme 
environmental conditions. We then described the six RT codes 
participating in this benchmark effort, focusing on how their im¬ 
plementation of the SHG calculation differs, presenting relevant 
heuristics for accelerating the calculation, and studying the ef¬ 
fects on the accuracy of the solutions. We also presented some 
practical hints with regards to the discretization of temperature 
and wavelength in the calculations. Most importantly, we pro¬ 
cessed the benchmark problem with each of the participating 
codes, and presented the results. 

We reported in detail on the similarities and differences be¬ 
tween the results from the participating codes and a reference so¬ 
lution. In the important wavelength range 3 jum < A < 1000/rm 
all participating codes reproduce the total dust emissivity within 
20% of the reference solution for all input fields used in this 
benchmark. Excluding the weakest and the softest input fields, 
the agreement in the same wavelength range is within 10%. 

Our discussion offered hints on how RT codes could be en¬ 
abled to properly calculate dust emission for a wider wavelength 
range. For example, when investigating systems with a lot of hot 
dust, such as circumstellar disks or accretion disks, it may be rel¬ 
evant to properly calculate dust emission for wavelengths shorter 
than 1 n m. To accomplish this, RT codes will need to model 
environment-dependent destruction of dust grains, and adjust the 
grain size distribution used in the dust emission calculation ac¬ 
cordingly. 

In conclusion, this benchmark effort shows that the relevant 
modules in RT codes can and do produce fairly consistent re¬ 
sults for the emissivity spectra of SHGs, which have a signifi¬ 
cant impact on the final result of a multi-wavelength RT simula¬ 
tion. We offer concrete, quantitative information on the level of 


(dis)agreement between RT codes, which will help inform the 
interpretation of RT simulation results that include SHG dust 
emission calculations of the type presented here. Specifically, 
this work paves the way for a more extensive benchmark effort 
focusing on the RT aspects of the various codes. And finally, 
we intend this work to serve as a reference for implementors of 
existing and new dust RT codes. 
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Fig. 1. The reference solutions generated with the public version of DustEM (see Sect. 5.1) using 3500 temperature bins and 250 iterations in the 
integral equation solver. The panels show the calculated dust emissivity for a selection of the input fields defined in Sect. 3. In each panel, the red 
curve represents the total emissivity and the other curves represent the portion of the emissivity for each grain type, silicate (magenta), graphite 
(green) and PAHs (blue). For the scaled Mathis input fields, the emissivity is divided by the input field strength U to allow identical axis ranges 
for all plots. 
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P. Camps et al.: Benchmark: Stochastic Heating of Dust Grains 



Fig. 2. A comparison of DustEM solutions for the most extreme input fields defined in Sect. 3, calculated with a varying number of temperature 
bins and iterations in the DustEM integral equation solver. The solutions employed as a reference for our benchmark are calculated with 3500 
temperature bins and 250 iterations; these solutions are represented in this figure by the zero-lines. The solutions calculated with the standard 
DustEM values of 200 temperature bins and 80 iterations are represented by the green curve. For these extreme fields, the standard solution 
deviates by up to 20% (and even more for wavelengths shorter than 1 pm). The solutions using 2500 temperature bins and 200 iterations (the blue 
curve) differ by less than 1% from the reference solution, indicating numerical convergence at these parameter values. The contribution of each 
grain type separately has a similar convergence behavior (not shown). 



A (pm) A (pm) A (pm) 

Fig. 3. A comparison of the emissivities calculated hy DustEM (using 3500 temperature bins and 250 iterations) for single-size, near-LTE grain 
populations to the corresponding equilibrium emissivities. The panels show the comparison for input fields ranging from extremely weak (left) to 
strong (right). The emissivity is divided by the input field strength U to allow identical axis ranges for all plots. We used a dust mixture consisting 
of 0.05 pm silicate (magenta) or graphite (green) grains with a total dust mass per hydrogen atom of 1(L 30 kg/H for each grain type. The solid 
curves represent the emissivities calculated by one of our codes (SKIRT) under the assumption of LTE. The dashed curves represent the solutions 
calculated by DustEM, without any LTE assumptions. The lower panels show the deviation of the equilibrium solutions from the corresponding 
full solutions. In a strong field, where we expect the grains to be in equilibrium, the solutions are indeed virtually identical. In a weaker field, the 
solutions differ since the grains are no longer completely in equilibrium. 
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A (|im) A (^m) 

Fig. 4. The relative differences between the total emissivities calculated by each of the codes participating in this benchmark and the corresponding 
reference solutions. The panels show the results for a selection of the input fields defined in Sect. 3. In each panel, the reference solution is 
represented by the zero-line. Positive percentages indicate results above the reference solution. The vertical dashed lines indicate where the 
reference solution becomes three orders of magnitude smaller than its peak value. 
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Fig. 5. The relative differences between the emissivities of the silicate component calculated by each of the codes participating in this benchmark 
and the corresponding reference solutions. The panels show the results for a selection of the input fields defined in Sect. 3. In each panel, the 
reference solution is represented by the zero-line. Positive percentages indicate results above the reference solution. The vertical dashed lines 
indicate where the reference solution becomes three orders of magnitude smaller than its peak value. 
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Fig. 6. The relative differences between the emissivities of the graphite component calculated by each of the codes participating in this benchmark 
and the corresponding reference solutions. The panels show the results for a selection of the input fields defined in Sect. 3. In each panel, the 
reference solution is represented by the zero-line. Positive percentages indicate results above the reference solution. The vertical dashed lines 
indicate where the reference solution becomes three orders of magnitude smaller than its peak value. 
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Fig. 7. The relative differences between the emissivities of the PAH component calculated by each of the codes participating in this benchmark and 
the corresponding reference solutions. The panels show the results for a selection of the input fields defined in Sect. 3. In each panel, the reference 
solution is represented by the zero-line. Positive percentages indicate results above the reference solution. The vertical dashed lines indicate where 
the reference solution becomes three orders of magnitude smaller than its peak value. 
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Fig. 8. The contributions calculated in the stochastic (solid) and equilibrium (dashed) regimes to the total emissivity (dotted), for each of the grain 
types silicate (magenta), graphite (green) and PAHs (blue), by one of our codes (DIRTY). The panels show the results for a selection of the input 
fields defined in Sect. 3. For the scaled Mathis input fields, the emissivity is divided by the input field strength U to allow identical axis ranges 
for all plots. The relative contributions for each regime may vary between codes because of the differences in the schemes to transition from one 
regime to the other. 
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Fig. 9. The grain size a^ns for which the participating codes transition from the stochastic to the equilibrium calculation regime. The panels show 
the smallest grain size that was considered to be in equilibrium for silicate (left), graphite (right) and PAH (right) grains, for each of the codes, 
across all of the input fields defined in Sect. 3 (within each panel, the scaled Mathis ISRF fields to the left, and the diluted black body fields to the 
right). 
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